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Abstract: We propose a new model for the description of complex granular 
particles and their interaction in molecular dynamics simulations of granular 
material in two dimensions. The grains are composed of triangles which are 
connected by deformable beams. Particles are allowed to be convex or concave. 
We present first results of simulations using this particle model. 



1. Introduction 



Molecular dynamics has been proven to be a well suited method to investigate 
the dynamic and static behavior of granular matter. The concept of molecular 
dynamics was initially used to simulate the dynamics of atoms and molecules, 
i.e. of particles without inner degrees of freedom. Alder and Wainwright who 
belong to the inventors of this method investigated already 1957 numerically 
the low density approximation of the entropy S = —N k B (In f(p,t)) (/ 

1 E-mail: thorsten@hlrsun.hlrz.kfa-juelich.de 

2 E-mail: volkhard@itp02.physik.hu-berlin.de 



Preprint submitted to Elsevier Science Publishers 



16 July 1995 



is the one-particle probability density in momentum space) in a hard sphere 
system consisting of N = 100 particles [1]. There are applications of molecular 
dynamics in many fields of interest and molecular dynamics became one of 
the standard methods in computational physics. Recent molecular dynamics 
simulations are of very high algorithmic complexity, and systems consisting of 
up to 10 9 particles (with simple interaction forces) have been simulated (see 
e.g. [2-4]). A very interesting overview on the historical development is given 
in the first part of Hoover's book [5]. 

In molecular dynamics of granular particles, sometimes called "granular dy- 
namics", each of the "molecules" is a macroscopic body typically with a di- 
ameter of the order D ps 0.1... 1 mm which has its own thermodynamic 
properties, i.e. it has internal degrees of freedom and hence it can dissipate 
mechanical energy when it is subjected to outer forces. In some special cases 
as in simulations of planetary rings [6,7] the diameter can also be of the order 
D rs 1 mm ...10m. In the most simple case one can express the inelastic na- 
ture of the collisions of the grains by introducing restitution coefficients which 
describe the relative velocities in normal and tangential direction v[- and 
of two granular particles i and j after a collision as functions of the relative 
velocities before the collision and V? 
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Here the relative velocity is given by 



(3) 



or, including the particle rotation with angular velocities Ui and u)j , by 

Vij = (?i + Ui x Ri (—nfj — (fj + ujj x Rjnj (4) 
with the unity vectors 
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Hence the relative velocities at the point of contact are 



Vij = n [Vij ■ n) 
vjj = tfy ■ t) . 
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In general the normal and tangential restitution coefficients and tx are 
themselves functions of the relative impact velocities = e N (v^ , V?^J and 

V; 1 ? , Vij)- F° r the case that the particles are rough homogeneous 



spheres the coefficients and e*r can be calculated analytically [8] solving 
the viscoelastic equations of the particles. Using this model one neglects the 
complicated interaction between particles during the collision. The approxi- 
mation is justified if collisions are relatively rare events, i.e. if the mean free 
path between collisions is long and one can restrict the calculation to two- 
particle-interactions only. If a system is in this regime one calls it a "granular 
gas". Event driven algorithms have been applied in many papers, recently e.g. 
in [9,10] for an investigation of inelastic clustering of granular particles (using 
constant restitution coefficients) and by Luding et al. to simulate vertically 
shaken material [11,12]. For sophisticated implementations of such algorithms 
see [13,14]. 



In many cases of interest, however, the particles collide very often or they touch 
each other permanently, i.e. the system reveals static properties. Examples for 
such systems are flows through pipes and hoppers, ball mills, shaken materials, 
sand heaps and many others. In the model proposed by Cundall and Strack [15] 
and Haff and Werner [16] the particles feel restoring forces in normal and 
tangential direction, and they lose kinetic energy during the collisions. If the 
particles i and j at the positions f and fj with radiuses Ri and Rj touch each 
other they feel the force 

F ij = Fgil + F$r (9) 

with 
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Y is the Young modulus, 7jv and 7t are the damping coefficients in normal 
and tangential direction and ji is the Coulomb friction coefficient. Eq. (12) 
describes the relative velocity of the surfaces of the particles at the point of 
contact and eq. (13) gives the effective mass. Eq. (11) takes the Coulomb 
friction law into account, saying that two particles slide on top of each other 
if the shear force overcomes fi times the normal force. Eq. (10) comes from 
the Hertz law [17] for the force between two rigid spheres which are in contact 
with each other. 



This model has been successfully applied to simulate the behavior of dry gran- 
ular material in many works describing various physical phenomena. It devel- 
oped to the standard method for calculations of granular dynamics. Some of 
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these problems are size segregation and convection in vibrated material in 
two and three dimensions (e.g. [18-21]), the flow in hoppers (e.g. [22-24]) 
and pipes (e.g. [25,26]), the flow in rotating cylinders (e.g. [27-29]), the mo- 
tion of granular material on an inclined surface [30,31], sound propagation in 
granular material [32], the onset of turbulence [33,34] and many others. Many 
experimental results of many authors could be reproduced numerically using 
this type of molecular dynamics simulations. There have been developed very 
efficient algorithms for the case of short range interaction with large force 
gradients as it is typical for granular materials, e.g. [35-37]. The model by 
Cundall and Strack yields reliable results in more dynamic systems, i.e. when 
the static friction of the particles does not play a major role. When the static 
properties of the material govern the behavior, however, the model might fail. 
In [38] is shown that one cannot build up a stable sand heap with a finite 
inclination of such particles but the heap dissolves under its own mass and 
the angle becomes smaller with rising number of particles. Lee [39] simulated 
static friction effects using spheres which are connected by springs. When two 
particles touch each other for a certain time a spring "grows" between these 
particles and keeps the particles together when forces are applied which would 
separate them. When the separating forces become too large the springs break 
and the particles move freely. Using this model Lee reproduced the finite angle 
of repose in a sand heap, however, the growing and breaking springs have no 
direct equivalent in nature, and hence the model seems to be quite artificial. 

Another simulation method which corresponds to the previous one for the case 
of very large damping has been introduced by Visscher and Bolsterli [40]. In 
each iteration step all of the grains are moved according to the motion of a 
vessel containing the granular material. Then the particles are released one by 
one, beginning with the lowest one (with respect to the direction of gravity). 
Each particle moves until it reaches the next local minimum of the potential 
energy when it hits either the wall or other previously released grains. When 
a particle reached its local minimum it stays in that position, i.e. a falling 
particle cannot cause the motion of another previously released particle. This 
behavior means that inertia of the moving particles is neglected as soon as they 
hit another particle or the wall, and hence it corresponds to the case of very 
large damping. The advantage of this algorithm is its numerical simplicity, 
i.e. it is possible to investigate much larger systems than using molecular 
dynamics. Recently it has been applied in simulations of several problems as 
size segregation [41] and the motion of particles in a rotating cylinder [42,43]. 
Obviously the simplification of infinite damping is not valid in each case, and 
hence one has to carefully investigate whether it is justified to apply this 
algorithm (see e.g. [45]). 

There are some more simulation methods for the dynamics of granular ma- 
terial which we want to mention only: Peng and Herrmann applied a Lattice 
Gas Automaton to investigate density waves in a pipe [46]. Caram and Hong 
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investigated granular flows using a Random Walk technique [47] . 



There are some models to simulate complex shaped grains which are composed 
of spheres. In the model by Gallas and Sokolowski [48] the grains consist of two 
spheres connected with each other by a stiff bar. Walton and Braun applied 
more complicated particles consisting of four or eight spheres rigidly connected 
with each other [49]. Using this model they examined the transition from sta- 
tionary to sliding flow and the transition from sliding to raining flow in a 
rotating drum. Poschel and Buchholtz [50,38] describe grains built up of five 
spheres where one of them is located in the center of the grain and four identi- 
cal spheres are at the corners of a square. Each pair of neighboring spheres is 
connected by a damped spring. The latter model was applied to the rotating 
cylinder [29] and it was shown that the simulation results agree much better 
with the experiment than equivalent calculations with spheres. Especially it 
was shown that one can reproduce stick-slip motion and avalanches which was 
not possible in simulations using spheres. The inclination of the surface of the 
material in the rotating cylinder and the dependency of the inclination on the 
angular velocity, however, did not agree well with the experiment. 

Mustoe and DePorter [51] proposed a particle model where the boundary 
geometry of the particles is defined (in local coordinates) by 
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Varying between 2 and oo the shape of the ith particles varies continuously 
from elliptical to rectangular. To detect contacts of pairs of particles described 
by eq. (14) one has to solve numerically equations of high order for each 
pair of possibly touching particles in each iteration step. This requires an 
iteration technique which converges slowly for higher exponents n« and which 
makes hence the algorithm numerically complicated. Hogue and Newland [52] 
investigated the flow of granular material on an inclined chute and through a 
hopper using a model where the boundaries of the convex particles are given 
by polygons with up to 24 vertices each. To detect whether two particles touch 
each other one has to calculate the intersections between each pair of vertices. 
During collisions energy is dissipated according to Stronge's energy dissipation 
hypothesis [53] in normal direction, whilst Coulomb's friction law models the 
energy losses also in tangential direction. 

Using grains which consist of interconnected spheres or of particles described 
by eq. (14) it is not possible to simulate particles with sharply formed corners. 
For some effects it seems to be essential to simulate such particles to reproduce 
the experimental observed effects. This point is discussed in detail in [38]. 
Tillemans and Herrmann [54,55] proposed a two dimensional particle model 
where the particles are convex polygons. When two particles i and j touch 
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each other, i.e. when there is an overlap, there acts the force 
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A is the compression area of the particles (overlap area), Y is the Young mod- 
ule of the material. The effective mass m e /^ is given by eq. (13), the relative 
particle velocity is v rel — fi — fj, where r*j points to the center of mass of the 
particle i, and L c is the characteristic size of the particle. The unit vectors are 
in the direction of the line that connects the intersection points of the overlap- 
ping particles i and j (e T ) and perpendicular to this line (e*). Only convex 
particles are allowed. This model was used e.g. in simulations of shear cells, 
earth quakes and flow through hoppers [55]. The results differ significantly 
from similar calculations using spheres. For the case of the hopper simulation 
they found clogging and arching which could not be found in simulations with 
spheres. A similar model was used earlier by Handley [56] who investigated the 
fracture behavior of brittle granular material. Potapov et al. [57,58] investi- 
gated a similar model for solid fracture. Initially they subdivide a macroscopic 
two dimensional body into many small equilateral triangles (elements). The 
forces in the body are resolved as forces at inter-element contacts of two dif- 
ferent types which they call "glued" and "collisional" . Glued contacts are the 
joints between elements interior to a solid body which can support stresses. 
When the tensile stresses exceed a certain value, cracks form and the glued 
bond breaks. Collisional contacts are contacts between the surface elements of 
the body (which eventually collides with other bodies or walls), and contacts 
between triangles in the bulk of the body where glued contacts have been bro- 
ken. In the three latter models [54-58] the convex grains have been Voronoi 
polyhedra in two dimensions [59]. 



In this paper we present a new particle model for molecular dynamics of 
granular material in two dimensions where the particles are simulated using a 
similar model which has some advantages. In the following section we describe 
the model, the forces acting between the grains are derived in sections 3 and 
4. Some aspects of the implementation and the performance of the algorithm 
are discussed in section 5 and some sample results based on this particle 
interaction model will be briefly discussed in section 6. 
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2. The model 



In our model the grains are composed of an arbitrary number of ideal elastic 
triangles which are connected by beams (see fig. 1). 




Fig. 1. Examples of grains composed of different numbers of triangles. The model 
is not restricted to convex grains. 

A beam in our sense is an deformable damped bar which is subjected to forces 
in the direction of its axis and transverse to its axis, i.e. to normal and shear 
forces, and to moments acting on its ends. Triangles which belong to the 
same grain do not interact with each other. When two triangles of different 
grains collide, i.e. if there is an overlap between both, they feel a restoring 
force. Hence the force acting upon the triangle % belonging to the grain j 
(je {1,2,. ..,7V}) is 

N n t (l) n J b (i) 

Fl= E E^+E^ + ^rt,^). (is) 

l=l,tyj k=l k=l 

The first sum in the first term in eq. (18) runs over all grains / except the 
jth, the second sum runs over all triangles k e {1, 2, . . . , n t (l)}, the /th grain 
consists of. rf k is the force which acts between the triangles i and k which 
belong to different grains j and / (see section 3). The sum in the second term 
runs over all beams connecting the triangle % of the grain j with other triangles 
of the same grain. A/ fc is the force induced by the beam k acting on the ith 
triangle of the grain j. They originate from the distortion of the beam that 
connects the triangle % with another one of the same grain. The acting forces 
and momenta will be discussed below in section 4. A similar model for the 
beam forces has been introduced earlier by Herrmann et al. for the fracture 
of disordered lattices [60]. The last term $(rf , r/) describes the action of an 
external force, as e.g. gravity, on the triangles. 

During its deformation a beam dissipates energy similar to a linearly damped 
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spring, i.e. proportional to its deformation rate (see section 4). To avoid time 
consuming calculations due to Steiner's law each beam is fixed at both ends 
in the center of mass points of the triangles it connects. When the beam is not 
deformed by eventually applied forces (the beam is at rest) it is perpendicular 
to the neighboring edges of the triangles which it connects. Fig. 1 shows some 
examples of grains. The only restriction concerning the number, the shape 
and the position of the triangles is caused by the condition that the beams 
are fixed in the center of mass. 

The next section 3 gives a detailed description of the force that acts when 
triangles collide, in section 4 the forces due to deformed beams are derived. 



3. The interaction of the triangles 

In this section we describe the calculation of the term f{ l k in eq. (18) originat- 
ing from the compression of two triangles % and k which belong to the grains 
j and I, respectively. When the triangles are compressed there exists a virtual 
"overlap area" which leads to an elastic restoring force T^. For simplicity here 
and in the following we drop the upper indexes of the variables. Areas of trian- 
gles XYZ and quadrangles WXYZ we denote by A(XYZ) and U(WXYZ). 
Fig. 2 shows the most frequently found type of a collision between the triangles 
AiBiCi and AkBkCk, i.e. a corner of one triangle deforms a side of another 
one. 

The absolute value of the restoring force is given by the shadowed area 
A (S^S^Afc) of the triangle Si^A^ times the Young module Y. Its direction is 
perpendicular to the intersection line S1S2 (Poisson hypothesis, see e.g. [23]). 
Hence resulting momenta acting upon the triangles AiBiCi and AkBkCk read 

M (AiBiCi) =HG l xf tk (19) 
M (A k BkCk) = HGk x fjw = -HGk x f ik . (20) 

The vector HGi points from the middle point of the line Si £2 to the center 
of mass of the triangle AiBiCi. 

Although the case shown in fig. 2 is the mostly occurring type of collisions 
there are several other types which treatment has to be discussed below. 

Classifying the possible interactions of the triangles we define five different 
types of collisions which are slightly different from those in [57]. The detection 
of overlapping polygons is a very important problem in computer graphics 
too, where one frequently has to decide whether two objects cover each other, 
and which of them is in the foreground or in the background, respectively. 
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Fig. 2. The "standard" type of a collision between two triangles AiBiCi and A^B^Ck. 
The forces ±I\fc are directed perpendicular to the intersection line S\S2- Moments 
act with respect to the center of mass points Gi and Gj. The absolute value of the 



interaction force 



r. 



ik 



is proportional to the shadowed intersection area A (SiS^-Afc) 



Therefore we found it helpful and inspiring first to have a look into advanced 
methods in computer graphics, e.g. [61,62]. In all figures 2-7 the overlapping 
area is drawn extremely exaggerated. Typically the overlapping area of two 
interacting triangles does not exceed a few tenth of a percent of the area of 
the triangles. 



A(S,S 2 A;) 




Fig. 3. The first type of collisions: There are two intersection points S\ and S2 ■ Both 
lie at the same edge of one of the triangles. 



(i) For the first type (fig. 3) there are two intersection points Si and S 2 lying 
at the same edge of one of the interacting triangles. One of the situations drawn 
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in fig. 3 corresponds to the standard collision (fig. 2). We calculate the area 
A (jSi^Aj) (shadowed area) given by the intersection points Si and S 2 and 
the included point A { . Because the overlapping area is definitely smaller than 
half the area of the triangles the force is proportional to the minimum between 
A (S 1 S 2 A i ) and A (A^d) - A (S 1 S 2 A i ) ' 



ik 



T ik = Y ■ min {A (S 1 S 2 A i ) , A (AB.Q) - A (S 1 S 2 A i )} . (21) 



The force acts perpendicular to the line between the intersection points SiS 2 . 

(ii) The second type (fig. 4) is characterised by two intersection points lying 
on different edges for both triangles. We calculate the areas A (S\S 2 Ai) (dark 
gray) and A (SiS 2 A k ) (light gray). The force is proportional to the overlap 



ik 



V lk = Y- [min {A (S 1 S 2 A i ) , A (A^Q) - A (S^AA} + 

min{A (S 1 S 2 A k ) , A (A k B k C k ) - A (S&A*)}] . (22) 



It acts perpendicular to the line between the intersection points as in the 
previous case. 




Fig. 4. The second type of collisions: There are two intersection points S± and S2 
lying on different edges for both triangles. 

(iii) The third type (fig. 5) is characterised by four intersection points, each 
pair of two points lies on one edge of both triangles. Therefore the third 
edge of both interacting triangles does not intersect the edges of the other 
triangle. There are two pairs of forces proportional to the quadrangular area 
□ (SiS^&t) of the overlap 
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(23) 
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one acting perpendicular to the line S1S2 and the other perpendicular to the 
line S1S4. The pre-factor 0.5 is due to two pairs of acting forces instead of one 
for the first two collision types. 




Fig. 5. The third type of collisions: There are four intersection points Si, S2, S3 and 
S4. Each pair of points lies on one edge of either the first or the second triangle. 

(iv) The fourth type (fig. 6) is characterised by four intersection points too, 
but here one of the triangles has intersection points at all its edges. Because 
interactions of this type occur extremely seldomly we did not implement an 
exact calculation for this case, but we calculate two interactions of type 1 with 
the intersection points S1S2 or S3-S4 respectively, instead of solving the cor- 
rect problem. In real simulation the fraction of this type occurs approximately 
10~ 5 of the number of collisions. Nevertheless one has to deal with these ex- 
tremely seldom events since otherwise they might cause problems from which 
the system cannot recover. If one does not care about these events usually 
the system gains in a few time steps after the event occurs a huge amount of 
kinetic energy, i.e. the system explodes. 




Fig. 6. The fourth type of collisions: There are four intersection points Si, S2, S3 
and S4. One of the triangles intersects with all its three edges. 
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(v) The fifth and last type of interaction (fig. 7) is characterised by six 
intersection points S\- S$. Equivalent to the forth type we do not calculate the 
exact interaction but substitute this calculation by solving three interactions 
of the first type with pairs of interaction points Si-S 2 , S^-S^ and S 5 -S 6 . 

According to the forces we derive forces in parallel to the axes of the 
Cartesian coordinate system and moments M acting with respect to the center 
of mass points Gj and Gk as described in equations (19, 20). 




Fig. 7. The fifth type of collisions: There are six intersection points S±,. . . ,Sq. 



4. Stresses in beams 



When torques and forces are applied to a beam the beam deforms. Since all 
deformations are assumed to be small compared with the size of the beam one 
can superpose the deformations originating from different applying forces and 
torques [63,64]. Fig. 8 shows an (infinitesimally) deformed beam with radius 
of curvature p. We find the approximation 

/ dO\ 1 , dx 
tan^-j= — -cfa or p=- (24) 

with tan(d©) = d<d. 

The dashed line ss is the neutral surface length of which does not change 
during the deformation. All other fibres below and above the neutral surface 
are in tension or in compression, respectively. The length of a fibre in distance 
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Fig. 8. An infinitesimally deformed beam with radius of curvature p. The fibres below 
or above the neutral fibre SS are in tension or in compression. For a beam with 
quadratic cross section b 2 the deformation leads to the resulting moment M = jt~ ■ 



y from the neutral surface is (p + y) dQ = (1 + y/p) dx and hence the strain 
reads e = y/p. With Hooke's law a = E • e one finds a = Ey/p. Finally the 
resulting momentum M is 



M 



r EI 

/ oydA = — (25) 

A P 



where the moment of inertia I — J y 2 dA depends on the geometrical shape of 
the cross section of the beam. In the present paper we assume that the beams 
have quadratic cross section of width b and hence / = 6 4 /12. When we express 
the radius of curvature p by 

1 d 2 v 

rs- (26) 

where v is the deflection of the beam from its initial position (with the ap- 
proximation tan© Ki 0). From eqs. (25) and (26) one finds the basic equation 
for the deflection of a beam: 

d 2 v ,, M 



This equation has to be solved to find the forces and moments which act 
when a beam is deformed. The integration constants are used to satisfy the 
boundary conditions. 

There are three different deformation modes for beams: elongation, shearing 
and bending. In the following three subsections 4.1, 4.2 and 4.3 we will discuss 
the deformation rules in detail. 
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4-1. Elongation of beams 



Since the exact notation of the acting forces and moments as vector functions 
of the coordinates of the triangles the beam connects, is not very instructive 
but confusing due to the length of the formulae expressions we discuss for 
simplicity of notation the deformation of the beams here and in the following 
in local coordinates. Instead of providing the exact vector notation we rather 
discuss absolute values. In all cases the directions in which the forces and 
moments act are very clear, moreover they are given for each case explicitly in 
the figures 9-11 and in figure A.l. For the implementation of the algorithm, 
however, one should note that it is necessary to transform the given expressions 
for the forces and momenta into the coordinate system of the triangles. These 
transformations require a considerable part of the computation time. 




Fig. 9. The shape and the location of a beam at rest. Its length is assumed to be 
L, its ends A and B lie on the x-axis. When the beam is not deformed there do 
not act any forces Fa, Fb or moments Ma, Mb- In the following figures 8-11 and 
in fig. 14 the x-axis is drawn as a thin horizontal line. 



Fig. 9 displays the shape and the location of a beam of length L at rest, 
i.e. when no deforming forces Fa, Fb or moments Ma, Mb at the ends A and 
B of the beam apply. The x-axis is drawn with a thin solid line, i.e. at rest 
the beam lies on the x-axis. 



For the elongation deformation we find according to the linear Hooke law the 
restoring force 

F el = \ AL\ ■ E. (28) 



4-2. Shearing 



We want to calculate the moments and forces Ma, Mb, Fa and Fb of a sheared 
beam as drawn in Fig. 10. For a beam of length L/2 which is fixed at one end 
(at x = 0) and where acts the shear force F (in negative y-direction) at the 
other end (at x = L/2) we find for the local moment M(x) = F ( |f — xj and 
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X 



(29) 



hence eq. (27) reads 

£7 V2 

with the boundary conditions v(0) = and v'(0) = 0. We find for the vertical 
deviation A/2 at the point x = L/2 



A/2 



Fir 



24 EI 

force which acts at the free end of the beam. 

|A 



(30) 




Fig. 10. When the beam in fig. 9 undergoes a shear deformation A there act the 
forces Fa and Fb and the moments Ma and Mb • The forces and moments are given 
in eqs. (22). 



One can assume that the sheared beam in Fig. 10 is built up of two of such 
beams of length L/2 as described above. Finally one finds 

• (3D 



12 EI 
and hence 



F A = -F B = F = ^^A (32) 
M A = M B = F ■ | = ^A . (33) 



4-3. Bending 



When a beam undergoes bending deformation (fig. 11), the resulting forces 
and moments can be found by superposing the forces and moments that would 
act if the beam would be bent at each side separately. Therefore we can restrict 
ourself on calculating the forces and moments according to the single sided 
deformation drawn in fig. 12. 
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(^A^ \r : ^ : ^t^ 

t L t 

Fa h" H F B 



Fig. 11. A bent beam. The forces Fa and Fb and the moments Ma and Ms can be 
calculated by superposing the forces and moments of two single sided bent beams 
(fig. 12), one of them bent at the side A, the other one at the side B. The results 
are given in eqs. (32) and (33). 

The deformation in fig. 12, however, can be understood as a deformation of a 
beam with free ends according to a single sided moment M (fig. 13) superposed 
by a moment M b (fig. 12) that assures that the angle Qb (fig- 12) results to 
zero. 




Fig. 12. A single sided bent beam. The forces and moments result from the super- 
position of the deformation drawn in fig. 13 and the restoring moment Mb (fig. 12) 
which assures the angle 0b to be zero. The equations for the forces and moments 
are given in eqs. (27) and (29). 



First we consider one sided bending of the beam (fig. 12). For a beam as drawn 
in fig. 13 where the moment M acts at one of its both sides one finds easily 



Since the beam in Fig. 12 has the angle Qb = we find the moments acting at 
the ends of the beam in fig. 12 by superposing the moment M B , which causes 
a virtual angle 6^*, fulfilling the condition + Q* B = 0: 

m b = ^(-ey = E. (36) 
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Fig. 13. If the ends of the beam are not fixed but freely moving, the applied moment 
M causes bending due to the angles Qa and Qb given in eqs. (23). 



The moment Mb causes the angle 



f L ML /o^ 

9 - = "6%7 = -12EI (37) 
and the resulting angle 6^ (fig. 12) is 

®A = e A + Q7 = ^- (38) 

Hence, finally we find for the different moments and forces for the bending 
deformation (fig. 12) 

AEI 

M A = -j—Qa (39) 
M B = -^©a (40) 

and with 

M A = -F B L - M B (41) 

F B = -F A = ^-®a ■ (42) 

(Note that the moments Ma and M B in eq. (41) have to be added as if they 
would apply to the same end of the beam.) 

If the beam is bent at both ends by angles 6^ and 6b (Fig. 11) the resulting 
moments Ma and M B and forces Fa and F B can be calculated by a super- 
position of the forces and moments according to two independent bending 
deformations of the type discussed in the previous case. 

The angle 6a generates the moments M\ and M B according to eqs. (39, 40) 

AEI 

M* A = — r Q A (43) 
2EI 

M B = ^Qa (44) 
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and the angle 6b causes the moments 



OTP J 

MX = ~—®b (45) 
AEI 

M** = -— e B • (46) 



Hence we find the resulting moments 

M A = M* A + M* A * = - (^©a + ^©b) (47) 
M B = M B + M* B * = - (^©a + ^©b ) , (48) 
and with M A + M B = —F B L the resulting forces 

F B = ~Fa = ^- (©a + e B ) • (49) 



Comparing the results from sections 4.2 and 4.3 one can show that each defor- 
mation of a beam can be expressed by bending and elongation, i.e. the shear 
deformation need not to be considered. The proof is given in Appendix A. 

The deformation of the beams is damped proportionally to the deformation 
rates 

M f = e A , (so) 

M B d) = 0b, (51) 

FA id) = -i(va-v b )-AB, (52) 
F^ (d) = 7 (va ~ vb) ■ AB, (53) 

where 7 is the damping coefficient of the beam material. As proven in the 
appendix the deformation of the beam is in linear approximation completely 
determined by the angles Q A and @b and the length L. Hence there is no 
damping force acting in shear direction. 



5. Implementation of the algorithm 



We implemented a molecular dynamics algorithm using the particle model 
described in the previous sections in FORTRAN. A Gear predictor corrector 
scheme of fourth order [65,66] was applied to integrate Newton's equations of 
motion. We have to proceed both the forces caused by the elastic deformation 
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of the triangles as well as the forces induced by the beams for each iteration 
step (eq. 18). Since every particle can interact with each other one causing 
a high algorithmic complexity, we applied a Verlet neighborhood list [67] to 
decrease the amount of computer time. As discussed below (see table 1) the 
calculation of the triangle intersections is the most time intensive part of the 
calculation. Hence we applied an advanced version of the neighborhood list 
method to keep the fraction of the computation of the triangle-interactions 
as small as possible. 

Our list method acts in two steps. In the first step we prepare neighborhood 
lists for each triangle to reduce the number of possible interactions which has 
to checked. Because we have to investigate all possible triangle interactions 
the time for this step rises as the square of the total number of triangles T ~ 
(N ■ n t ) 2 . (For simplicity we assumed for this estimate that all N grains consist 
of n t triangles each. Hence there are (N ■ n t ) ■ ((N — 1) • n t ) possible triangle- 
interactions.) In the second step we check for each entry of the constructed 
neighborhood list whether the corresponding triangles interact indeed or not. 
If they interact, i.e. if there is an overlap area, the type of interaction is 
determined according to the classification scheme in section 3. 

Finally we construct three lists for the interactions of the first three types 
of the classification scheme. The interactions of the forth and fifth type are 
replaced by interactions of the first type as described above. The computer 
time for this step rises linearly with the number of particles N and with the 
number of triangles n t belonging to each particle. Different from usual Verlet 
tables these lists contain only triangles which definitely do interact, and hence 
we call the most time consuming subroutines for the calculation of the forces 
only for triangles which do interact in one of the five manners described above 
(section 3). 

The calculation of the forces induced by the beam deformations follows straight 
forward the formulas described in section 4. 

In table 1 we present a detailed performance analysis of the algorithm. The 
numbers in the table refer to simulations on a DEC 3000/700 workstation. As 
visible from the data most of the time (approximately 70%) is necessary to 
construct the interaction list, because we check here every pair of neighbored 
triangles for intersections. The construction of this lists give the possibility to 
classify the interactions and therefore to use for each kind of interaction its own 
optimised algorithm. This prevents us from calculating useless intersection 
points, areas and other data. 

The computer time needed for the predictor-corrector integration is very low 
(3.7%). Thus it seems to be not useful to optimise the integration procedure, 
or to decrease the accuracy of the integration to safe computer time. The 
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calculation of the beam forces is not time expensive as well. 
Table 1 

The performance analysis of the algorithm. The calculation has been done on a 
DEC-3000/700 workstation (explanation in the text). 



Basic algorithm 


portion of 
computation 


subroutine 


portion of 
computation 


Construction of the lists 


73.7 % 


Verlet list 


2.07 % 






Classifying interaction 


70.62 % 


Calculation of the 


7.87 % 


TYPE 1 


7.54 % 


triangle interaction 




TYPE 2 


0.24 % 






TYPE 3 


0.09 % 


Calculation of the 
beam forces 






13.33 % 


integration 


3.71 % 


predictor 


1.16 % 






corrector 


2.55 % 


all others 






1.4 % 



From the data in table 1 it is obvious that triangle interactions of the first 
type dominate. Approximately 98% of the occurring interactions are of this 
type. Interactions of the fourth and fifth type, where our program yields ap- 
proximative results only, occur extremely seldomly. We found them only a few 
times during all our simulations. 



6. First results 

6.1. The collision of a two particles 

To demonstrate the non-trivial behavior of colliding non-spherical particles 
and to check the correctness of our implemented model we want to present the 
results of an experiment where a quadratic grain collides with another resting 
one in the absence of gravity. The grains consist of four equal triangles each, 
forming a simple quadratic grain as shown in fig. 1. The parameters of the 
materials are: Y — 2 ■ 10 7 g/(cm sec 2 ), E = 1 • 10 5 g/ {sec 2 ), 1=1- 10~ 4 cm 3 
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and 7 = 9 gj sec. The value of the Young-module Y agrees with most of the 
simulations of spheres which can be found in the literature. 



Figs. 14 (left figures) show the time evolution of the particles as a sequence of 
stroboscopic snapshots for different relative velocities of the particles a) v re i = 
10 cm/ sec, b) v re i = 50 cm/ sec and c) v re i = 100 cm/ sec. 

One of the triangles of each grain has been filled to visualise the rotation of 
the grains. Initially the velocity of the left particle is v = (v re i, 0) cm /sec and 
the right particle rests. For low velocities v < v tc ~ 70cm/ sec the grains col- 
lide twice within a very short time. Therefore the traces of the particles that 
undergo soft collisions with velocity values v < v tc differ qualitatively from the 
traces of particles colliding with higher relative velocity. Since it is hard to dis- 
play the complicated motion of the grains in detail we refer to our World Wide 
Web-site URL http://summa.physik.hu-berlin.de:80/~thorsten/MDASGP.html [68] 
where one can find animated sequences of various problems mentioned in the 
current paper. The right hand side in fig. 14 shows the rotational and kinetic 
energies and the total energy as a function of time. The initial kinetic energy 
of the translational motion of the left particle results in kinetic energy of both 
particles due to translational motion and rotation. A part of the mechanical 
energy is lost due to dissipation in the beams. The relative amount of rota- 
tional energy depends strongly on the impact velocity as expected from the 
discussion above. The dependence of the extinction of the rotational degree of 
freedom on the impact rate is demonstrated in fig. 15 too. This figure shows the 

time evolution of the angular velocity of the grains u> = — 



nt 



nt 

E {Vi - Vg) X Sj 

i=l 



- 1 nt -> 

over time for different values of the impact velocity. v g — — E i>j is the ve- 

locity of the grain and Sj is the vector pointing from the center of mass point 
of the grain to the center of mass point of the triangle i. n t is the number of 
triangles the grain consists of. 

For high velocity one observes an abrupt change of the rotational motion at the 
time of the collision and a short damped oscillation according to the excited 
inner degrees of freedom of the grain (distortion of beams) and their relax- 
ation. After the collision the particles move on with constant translational 
and rotational velocities. For low impact rate we find a quite different behav- 
ior: First the particles collide as in the previous case, resulting in a certain 
rotational motion. Short time after the first collision, however, they collide a 
second time with different edges. Finally we find for the latter case of slow 
collisions a very small resulting angular velocity because the second collision 
causes angular moments in opposite directions, hence the resulting moment is 
small. 



Observing the motion of a moving particle colliding with a resting one (see 
i.e. the simplest possible contact of particles, one can remark that even for 
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Fig. 14. The collision of a horizontally from left to right moving grain (black 
drawn) with a resting one (grey) for different impact rates v re i = 10 cm/ sec (top), 
v re i = 50 cm I sec (middle), and v re i = 100 cm /sec (bottom). The stroboscopic snap- 
shots of the grains are take at equidistant time intervals. An animated sequence 
of the collisions can be accessed via World Wide Web [67]. The right hand figures 
display the kinetic, rotational and total energies of the particle collisions over time. 
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Fig. 15. The angular velocities of the particles over the time. For lower velocities 
the particles collide two times, therefore those collisions results in a small angular 
velocity of the grains. 

this system there is a complex behavior. For rigid spheres the velocities after 
a collision as a function of the impact velocities (eq. (1, 2)) can be calculated 
analytically solving the viscoelastic equations of the spheres [8] . We guess that 
this will not be possible for the case of our particles. 

6.2. Outflow of a hopper 

The outflow of a hopper is of high interest not only because of the technological 
importance of this process, but also because of the exciting phenomena as 
clogging and density waves which have been discovered recently (e.g. [69]). 
Density waves and clogging have been investigated using molecular dynamics 
in two and three dimensions in various papers, e.g. [23,70,71]. There is at 
least one other interesting and surprising phenomenon, detected recently by 
Evesque and Meftah [72,73]: They found that a hour glass "ticks" much slower 
if it is subjected to vertical vibrations y = A cos (cut). The effect seems to 
depend on the frequency / = ujjl-n and on the acceleration of the vibration 
T = Auj 2 . For frequencies between 40 Hz < f < 60 Hz surprisingly they 
observed for some accelerations that the flow almost stops. 

We investigated this effect using different particle models and we found that 
the effect could be reproduced with the new particle model only. The results 
will be described in detail in a forthcoming paper [74]. 

Fig. (16) shows snapshots of the system with N = 400 complex particles. For 
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an animated sequence we refer to our WWW-site [68]. The left figure shows 
the outflowing hopper, the flow varies irregularly and the grey scale codes for 
the particle velocities. When animating the figures one clearly observes density 
waves [68] . The right hand figure shows the same system a few moments later. 
The flow has stopped due to clogging. 




r.. 

Fig. 16. Snapshots of the simulation of an outflowing hopper. The left figure shows 
regular outflow, the right figure show a clogging situation. An animated sequence 
is available via World Wide Web [67]. 



6.3. Granular flow in a rotating cylinder 

The flow of granular material in a rotating cylinder is one of the most popular 
problems in the field of granular materials. It has been investigated experimen- 
tally and theoretically by many authors using various techniques (see e.g. [75- 
78,27,43,50] and many others). The results shall not be discussed here. We 
applied the algorithm to this problem too. For the detailed description of the 
results see [79]. Here we only want to demonstrate the abilities of the method 
and to mention the main results briefly. 

Fig. 17 shows snapshots of the simulation of a slowly rotating cylinder. The 
grey scale codes for the velocity of the grains, black means high velocity, 
black codes for high velocity. For low angular velocity of the driven cylinder 
one finds experimentally stick-slip flow (e.g. [75]), i.e. the material moves 
downwards not homogeneously but it forms avalanches. In the right figure 
one observes an avalanche on the top of the material indexed by light grey 
shadowed grains. The left figure shows the same system immediately before 
the avalanche. Avalanches are relatively seldom events. Most of the time the 
particles rest with respect to each other, i.e. they move only due to the rotation 
of the cylinder (left side of fig. 17). When increasing the angular velocity of 
the cylinder there is a relatively sharp transition between the stick-slip flow 
and the homogeneous regime. This transition could be found in the simulation 
using non-spherical grains [79]. Animated sequences for both regimes, stick- 
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slip and continuous flow, can be accessed via World Wide Web [68]. 




Fig. 17. Snapshots of the simulation with N = 500 complex particles with 
P £ [0.1 cm; 0.2 cm] in a rotating cylinder of diameter D = Acm. The angular 
velocity is $7 = 0.1 sec -1 . The wall particles have not been drawn. The right snap- 
shot shows an avalanche, the left one has been taken short time before. The grey 
scale codes for the particle velocity. 



Fig. 18 shows the dependence of the difference between the inclination of 
the material and the angle of repose on the angular velocity of the cylinder. 
The line displays the behavior found by Rajchenbach — @ c ~ Q 2 [75]. The 
simulation data are in good agreement with this observation. 
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Fig. 18. The inclination of the material surface over the angular velocity Q. 
The dotted line displays the function which has been measured experimentally by 
Rajchenbach [73]. 



Both results mentioned, the sharp transition between the flow regimes and 
the inclination as a function of the angular velocity of the cylinder, agree with 
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experimental observations. We have not been able to find these effects using 
simple spheres or using particles composed of spheres [50,29] in molecular 
dynamics simulations. 



7. Conclusion 

We presented a model for the simulation of the particles of a granular ma- 
terial. Each particle consists of triangles which are connected by deformable 
beams. There are no restrictions concerning the shape of the grains (convex or 
concave), nor the number or the shape of the triangles a grain consists of. To 
preserve calculation time it is favourable to chose arrangements of the trian- 
gles so that the beams are fixed at the center of mass points of the triangles. 
Otherwise one needs additional computation time caused by the Steiner law. 

When two triangles of different grains collide there acts an elastic restoring 
force proportionally to the compression of the particles. During the collision of 
the triangles the energy is preserved. A collision of a triangle with another one 
causes moments and forces acting on this triangle and hence a deformation of 
the beams which connect the triangles to neighboring ones. The deformation 
of the beams causes forces and moments acting on the neighboring triangles. 

Beams can be deformed in axial direction (elongation) and in shear direc- 
tion and they can be bent. When a beam is deformed it dissipates energy 
proportionally to the deformation rate. 

In the present paper we described beams that recover completely while dis- 
sipating mechanical energy when the deforming moments and forces vanish. 
Our model, however, is not restricted to this case. Interesting phenomena are 
plastic deformation and wear which can be easily simulated by introducing 
thresholds for the beam forces and moments or for the deformations. When 
a force or a moment or a deformation, respectively, exceeds the threshold the 
beam breaks or deforms permanently. Similar simulations with other models 
have been reported in the literature [55-57]. 

The proposed algorithm has been implemented in FORTRAN. We have demon- 
strated the behavior of a system of grains applying it to three examples of 
granular assemblies. In a simple collision simulation of two grains we dis- 
cussed the trajectories of the grains as well as the evolution of the kinetic and 
rotational energies depending on the impact rate. Sample results for the flow 
out of a hopper and the motion of granular material in a rotating cylinder 
have been presented. In the latter case we found quantitative agreement of 
the simulation with the experiment that could not be found so far, neither 
with spherical grains nor with non-spherical grains of other type. 
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Appendix A. Shear deformation can be expressed by bending 

The deformation of a beam in the two dimensional space has three degrees 
of freedom: the distance of the relative positions of the ends of the beams 



L — y {xa — xb) + {va — Ub) , and the angles Oa and 6#- These values are 
the parameters in eqs. (47, 48) and (49). The shear parameter A in eqs. (32, 
33) is a fourth parameter, i.e. one of them must be redundant. In the following 
we will show that shear deformation can be expressed by bending, at least in 
our approach of a linear deformation-force relation. 



Fig. A.l. The deformation of a beam in shear direction A (upper part) can be repre- 
sented by bending deformation (lower part). The appendix contains the proof that 
one need not to consider shear and bending deformation in the linear approximation 
which is justified for small deformations. 

Suppose we have a deformed beam as drawn in Fig. A.l (bottom). Then 
applying eqs. (32, 33) we find for the forces and moments according to shearing 






(A.l) 



n 




12EI 



(A.2) 
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and according to bending (eqs. (47,48) and (49)) 



9FJ 

M\ = -—Q* B (A.3) 

M b B = -^Q B (A.4) 

F b A = -F b B = -^^s- (A.5) 
The total forces and moments read 

M\ = M% + M\ = ^ A - ^ 0^ (A.6) 

M* = M B + M b B = ^ A - ^ e% (A.7) 

(L*) (L*) B V ; 

Turning the entire system by the angle Qa ~ tanO^ = — we obtain the 
new system as drawn in the top of fig. A.l with the new angles Qa = — j; and 
Ob = Q* B — j;. With eqs. (47, 48) and (49) we find 

M, = -^e,-^e B = ^A-^e^ (a.9) 

L 



2EI 


e B = 


6EI 




2EI 




A- 




L 




L 2 




L 


AEI 


e B = 


QEI 




AEI 




A- 




L 




L 2 




L 



M B = -^e,-^pe B = ^A-^0^ (a. io) 

FA = - FB = - 6 ^L ( e A + e B)= 1 JELA- 6 -§Le' B (a.h) 

Approximating L by \JL 2 + A 2 = L* evidently the results in eqs. (A.9-A.11) 
coincide with eqs. (A.6-A.8). Thus we do not need to consider shear deforma- 
tion in our simulations since it is redundant here. 
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